Rainfall vs Niño#

https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/doc/GHCND_documentation.pdf

import warnings
warnings.filterwarnings("ignore")
import os
import os.path as op
import sys

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import df2img


sys.path.append("../../../../indicators_setup")
from ind_setup.colors import get_df_col, plotting_style
from ind_setup.tables import plot_df_table
from ind_setup.plotting_int import plot_oni_index_th
from ind_setup.plotting import plot_bar_probs_ONI

plotting_style()
from ind_setup.core import fontsize

sys.path.append("../../../functions")
from data_downloaders import GHCN, download_oni_index

Define location and variables of interest#

country = 'Palau'
vars_interest = ['PRCP']

Get Data#

update_data = False
path_data = "../../../data"
if update_data:
    df_country = GHCN.get_country_code(country)
    print(f'The GHCN code for {country} is {df_country["Code"].values[0]}')

    df_stations = GHCN.download_stations_info()
    df_country_stations = df_stations[df_stations['ID'].str.startswith(df_country.Code.values[0])]
    print(f'There are {df_country_stations.shape[0]} stations in {country}')

Using Koror Station#

if update_data:
    GHCND_dir = 'https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/access/'
    id = 'PSW00040309' # Koror Station
    dict_prcp = GHCN.extract_dict_data_var(GHCND_dir, 'PRCP', df_country_stations.loc[df_country_stations['ID'] == id])[0]
    data = dict_prcp[0]['data']#.dropna()
    data.to_pickle(op.join(path_data, 'GHCN_precipitation.pkl'))
else:
    data = pd.read_pickle(op.join(path_data, 'GHCN_precipitation.pkl'))

st_data = data

ONI index#

https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php

p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
if update_data:
    df1 = download_oni_index(p_data)
    df1.to_pickle(op.join(path_data, 'oni_index.pkl'))
else:
    df1 = pd.read_pickle(op.join(path_data, 'oni_index.pkl'))
lims = [-.5, .5]
plot_oni_index_th(df1, lims = lims)
st_data_monthly = st_data.resample('M').mean()
st_data_monthly.index = pd.DatetimeIndex(st_data_monthly.index).to_period('M').to_timestamp() + pd.offsets.MonthBegin(1)
st_data_monthly
PRCP
DATE
1947-08-01 NaN
1947-09-01 NaN
1947-10-01 NaN
1947-11-01 NaN
1947-12-01 NaN
... ...
2024-10-01 13.756667
2024-11-01 7.055172
2024-12-01 7.751852
2025-01-01 17.175000
2025-02-01 8.452632

931 rows × 1 columns

rolling_mean = 6 #months
df1['prcp'] = st_data_monthly['PRCP'].rolling(window=rolling_mean).mean()
fig, ax = plt.subplots(figsize = [15, 6])
df1.ONI.plot(ax = ax, color = get_df_col()[0], lw = 2)

ax2 = ax.twinx()
df1.prcp.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tmax.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tdiff.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tmean.plot(ax = ax2, color = get_df_col()[1], lw = 2)
<Axes: >
../../../_images/7b63303e8f988f717156269cca1076387ef5d95a83f9e37a66079cc94ea93672.png
low_lim = np.nanmin(df1.prcp)

fig, ax = plt.subplots(figsize = [15, 6])
df1.prcp.plot(ax = ax, color = get_df_col()[1], lw = 2)

ax.fill_between(df1.index, low_lim, df1.prcp, where = (df1.ONI > lims[1]), color = 'lightcoral',
                 alpha = 0.7, label = f'ONI over th: {lims[1]}')
ax.fill_between(df1.index, low_lim, df1.prcp, where = (df1.ONI < lims[0]), color = 'lightblue', 
                alpha = 0.7, label = f'ONI below th: {lims[0]}')

ax.fill_between(df1.index, low_lim, df1.prcp, where = ((df1.ONI > lims[0]) & (df1.ONI < lims[1])),
                 color = 'lightgrey', alpha = 0.075)

ax.legend(fontsize=fontsize)
ax.set_title('PRECIPITATION and ONI', fontsize = fontsize)
ax.set_ylabel('PRECIP', fontsize = fontsize)
ax.set_xlabel('Time', fontsize = fontsize)
Text(0.5, 0, 'Time')
../../../_images/85010efd042a34f3a6b15f4b160f2350b85121f6895067d110f1e4bfb040bfd9.png
plt.figure(figsize=(5, 4))
sns.heatmap(df1.corr(), annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('Correlation Heatmap')
plt.show()
../../../_images/13e9e1460aa24182892d3d383ecc8378df5cbd2f83b062b3f1ed712a3b84e870.png
df2 = df1.resample('Y').mean()
df2['oni_cat'] = np.where(df2.ONI > lims[1], 1, np.where(df2.ONI < lims[0], -1, 0))   
fig = plot_bar_probs_ONI(df2, var='prcp', y_label = 'Mean Annual Precipitation')
../../../_images/9474e5ba29f95d6f5958c43c1c55fbd97afd559958dc91adcb66defcbd233bca.png
df2['prcp_ref'] = df2.prcp - df2.loc['1961':'1990'].prcp.mean()
fig = plot_bar_probs_ONI(df2, var='prcp_ref')
fig.suptitle('Precipitation Anomaly over the 1961-1990 mean', fontsize = fontsize)
Text(0.5, 0.98, 'Precipitation Anomaly over the 1961-1990 mean')
../../../_images/c977b25d71be3b51f5b0b5a9e066210f7765ed1ac4c0f7e47dfd869a93cb7d16.png
df_format = np.round(df1.describe(), 2)
fig = plot_df_table(df_format, figsize = (400, 300))
df2img.save_dataframe(fig=fig, filename="getting_started.png")
../../../_images/f8241bb32104758f7a6251d879b45f6a140bbd156fc345691c8f9d3dceaca05e.png